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We demonstrate that the phase response curve (PRC) can be reconstructed using a weighted 
spike-triggered average of an injected fluctuating input. The key idea is to choose the weight to 
be proportional to the magnitude of the fluctuation of the oscillatory period. Particularly, when a 
neuron exhibits random switching behavior between two bursting modes, two corresponding PRCs 
can be simultaneously reconstructed, even from the data of a single trial. This method offers an 
efficient alternative to the experimental investigation of oscillatory systems, without the need for 
detailed modeling. 
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Synchronization phenomena in networks consisting 
of interacting nonlinear dynamical elements exhibiting 
limit-cycle oscillations have been the subject of inten- 
sive study in physical, biological and social systems 
In general, such a system is described by ordinary dif- 
ferential equations of the form — F(X), where X 
denotes a multi-dimensional state of the system. Ow- 
ing to nonlinearity and complicated interactions, how- 
ever, we often encounter difficulties when attempting 
to analyze synchronization properties of such systems 
(see arrow A in FigQ]). One of the most successful ap- 
proaches for dealing with such difficulties is the phase 
reduction method (arrow B), in which the evolution of 
each active element is described by only one degree of 
freedom, the phase [2]. With this method, the dy- 
namics of N coupled oscillatory systems can generally 
be described by a set of equations of the form = 

u>i + X)j=i Tij (4>i — (j)j), where u>i is the natural frequency 
without coupling. The phase of the z-th system repre- 
senting the timing of its limit-cycle oscillation is denoted 
by 4>i . This description is valid if the perturbation due to 
the interactions is so weak that its effect only changes the 
phase asymptotically. The coupling function Tij(<j>) can 
then be derived from the original dynamical system as 
Z l {d)V l3 {X^ ) {e),X^{9-mO, where 
Vij{Xi, Xj) and Xq\<P) represent the interaction ex- 
erted by the j-th system on the i-th system and the limit- 
cycle solution for the uncoupled i-th system, respectively 
[3(. Here, Zi represents the phase response curve (PRC) 
of the i-th. system. Therefore, if we can obtain the PRC, 
then we can predict the synchronization properties of the 
coupled system for any weak interaction (arrows C). 

It is well known that the PRC can be calculated from 
the original model equations with the adjoint method 
Owing to technical difficulties limiting the applicability 
of experimental studies, however, it is often difficult to 
construct system-specific detailed models including the 
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FIG. 1: Without detailed modeling, the dynamical system 
reconstructed from the PRC obtained through direct exper- 
imental measurements can be used to correctly predict the 
collective synchronization properties of a set of coupled sys- 
tems. 



essential properties for synchronization (arrow D). For 
instance, if we wish to determine how a given neuromod- 
ulator causes the synchronization properties of a coupled 
neuronal system to change, we cannot evaluate all the 
effects of the neuromodulator on the dynamics of the 
various ion channels. As an alternative approach, there 
have been some recent attempts to directly measure the 
PRC experimentally (arrow E). The most important ad- 
vantage of this approach is that, once the exact PRC is 
known, the dynamics of a set of weakly coupled oscil- 
lators are fully determined, whether or not the original 
detailed dynamics are known. In the example considered 
above, therefore, all the effects of the neuromodulator 
are reflected in the measured PRC, and with this, the 
dynamical behaviors can be correctly predicted. 

The standard method for obtaining PRCs through di- 
rect experimental measurements is to measure the phase 
shift induced by stimulating the system with a brief 
pulse-like perturbation at various phases of the period 
5]. Unfortunately, in such a study we often face the dif- 
ficulty that the information from which we can derive 
the PRC is lost in the noise arising from the uncontrol- 
lable nature of the environment, though there are several 
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FIG. 2: Schematic of the experimental procedure to recon- 
struct the PRC in our method. (a) During the experi- 
ment, the injected fluctuating current I(t) and the ISIs n 
are recorded. The recorded I(t) is divided into the current 
traces hit) in the individual intervals, (b) The rescaled cur- 
rents Ii(ti) are obtained by adjusting the duration of each 
current trace U(t) to the average T. (c) In the WSTA, the 
PRC can be reconstructed by averaging over all the rescaled 

'<■, the normalized 



currents Zj(ij) with the weights Aj = - 
difference between the ISI and the average interval. 



statistical fitting methods which have been applied to ef- 
ficiently extract the PRC from noisy data [6|. 

An alternative approach for obtaining the PRC which 
we consider here, is to design a different experimental 
procedure for the measurement, instead of the standard 
method described above. In connection to this, G. B. 
Ermentrout et al. have recently shown that the spike- 
triggered average (STA) of the injected noisy current is 
proportional to the derivative of the PRC under suitable 
conditions Q- This result inspired us to develop a more 
useful, novel method for experimental protocols. 

In this paper, we demonstrate that the PRC can be 
directly reconstructed using an appropriately weighted 
spike-triggered average of the injected fluctuating inputs. 
The key idea in this method is to choose the weight to be 
proportional to the magnitude of the fluctuation of the 
oscillatory period. In the conventional method employing 
pulse-like perturbations, the spontaneous fluctuation of 
the period tends to make the estimation of the PRC quite 
rough. Interestingly, the disadvantage of the fluctuating 
period in the case of the conventional method becomes 
an advantage in our method. Although the proposed 
method can be applied to various real oscillatory sys- 
tems, such as Belousov-Zhabotinsky reactions, for sim- 
plicity, we consider two neuronal systems to explain the 
experimental procedure in the following. 

Experimental procedure - We consider the situation in 
which a neuron with some constant bias current exhibits 
spikes regularly with a period T. When an additional 
fluctuating current with zero mean is injected into this 
neuron, the inter-spike interval (ISI) generally fluctuates 
about the average value T, as shown in Fig^a). As the 
first step of the experimental procedure, we record the 
stimulus current I(t) and all the spike times over the en- 
tire spike train. Next, we divide the stimulus I(t) into in- 
dividual stimuli Ii(t) between successive spikes, whose in- 
terval duration is denoted by t$ (i = 1, . . . , N; N denotes 
the number of the recorded spikes). Then (FigJ^b)), the 
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FIG. 3: The estimated PRCs obtained from the WSTA in the 
case of the Hodgkin-Huxley model, (a) The time evolution of 
the membrane voltage (upper solid trace) and the injected 
fluctuating current (lower trace). The injected current con- 
sists of two components: the bias current 7bi as = 10^A/cm 2 , 
and the zero-mean O-U noise with a — 0.5/iA/cm 2 . In the 
case that there exists only bias current, the neuron exhibits 
regularly periodic firing with a frequency of approximately 
68Hz (upper dotted trace), (b) The WSTA-estimated PRCs 
for three different measurement durations, Ltrial- The esti- 
mated PRC rapidly converges to the correct one as Ltrial in- 
creases. 



slightly different duration times of the stimulus currents 
Ii(t) are rescaled to the uniform period T as 



U 



T 



-t. 



U e [0,T] 



(1) 



Finally, the weighted spike-triggered average (WSTA) of 
the rescaled stimulus current Ii(ti) is defined as 



WSTA(fi) = (Ailifa)) 



(2) 



where Aj = ~ Ti , and the angular brackets represent 
the average over all spikes, i.e. the index i (FigJ^c)). 
As shown below, we can prove that the above-defined 
WSTA is proportional to the PRC, provided that the 
fluctuations in the stimulus current are not too strong. 

Theory - When a neuron fires almost periodically under 
the influence of a zero-mean fluctuating stimulus current 
I(t), the corresponding phase dynamics can generally be 
described by the equation 



(3) 



in which u) = and Z(<f>) represents the PRC with re- 
spect to the input current. Let us consider the effect of 
the z-th stimulus, Ii{t) = (j£i(t), where the parameter 
/x is introduced to account for the effect of the magni- 
tude of the fluctuating current. Hereafter, we omit the 
subscript i of the rescaled time U for simplicity. In ad- 
dition, we assume that is a stationary stochastic 
process with unit strength. The auto-correlation of 
is then defined as C(t) = (£i(?)£i(0)). Then, substituting 
I(t) = fi£i(t) into Eq.([3]), and integrating from t = to 
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Tj, we obtain Aj = ^ J Q T Z(</) i (s))^(s)ds, where Eq.fT]) 
has been used. Next, substituting this into Eq.([2]), we 
find 

WSTA(t) = £- / T <Z(^(5))6(5)6(t)) d5. (4) 
27r Jo 

We assume that we can expand the time development of 
the phase with respect to fi as (f>i(t) = cot + ficf^it) + 

\J?<\>f" > (£) + •" ■ Using this and the Taylor expansion in 
Eq.fll]), we get 

WSTA(i) 



u 2 r T 

£- / Z(w5)C(s-i)dS + 0(/x 3 ). 
2tt Jo 



(5) 



Assuming that the timescale on which the auto- 
correlation of the stimulus current decays is much shorter 
than the period T, we can approximate C(t) as 8(t), 
where 8 is the Dirac delta function. Using the fact that 
fi is then given by (li(t)li(0)) = fi 2 8(t), we finally obtain 



WSTA(i) = £-Z(^) + 0(// 3 ), 

Z7T 



(6) 



where <fi = tot. We have thus found that the PRC can 
be reconstructed from the spike-triggered average of the 
stimulus current using a weight proportional to the nor- 
malized difference between the ISI and the average inter- 
val. 

In the numerical examples considered below, to gener- 
ate fluctuating stimuli, we used the zero-mean Ornstein- 
Uhlenbeck (O-U) process dl t = --fldt + y/2-fa 2 dW t , 
where Wt is a Wiener process. For sufficiently large 7 
and in the long-time limit, we can approximate (Itlo) 



as cr 2 e 7 '*' 



a: 



a 2 e 



J\t\ 



dt)8(t) 



-8(t) (i.e. \l 



y/2/j a). In all the examples considered here, we fixed 
7 = 5 and varied a as the control parameter. 

Numerical examples - To confirm the validity of our 
proposed method, we compare the PRCs obtained from 
the WSTAs and the PRCs calculated using the adjoint 
method. In FigJ31 we first show typical results in the case 
of the Hodgkin-Huxley model |8j. Figure a) displays 
a typical voltage trajectory in response to the fluctuat- 
ing current. Figure [3fb) presents comparisons between 
the true PRC (black trace) and the PRC estimated from 
the WSTA (red trace) for three different measurement 
durations. We see that as the measurement duration in- 
creases, the WSTA-estimated PRC converges to the true 
one. This result suggests that the PRCs of real systems 
can be accurately estimated within practically reasonable 
recording times experimentally. 

We next consider the more complicated situation in 
which a bursting neuron exhibits singlet or doublet firing 
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FIG. 4: From a single fluctuating current (Ltriai =180 sec) 
that causes a bursting neuron to exhibit irregular singlet 
or doublet firing, two PRCs corresponding to the singlet 
and doublet states can be reconstructed simultaneously, (a) 
The membrane voltage of the bursting neuron (upper trace) 
and the injected current (lower trace, a = 0.2/^A/cm 2 and 
-^bias = 1.605/iA/cm 2 ). (b) Distributions of conditional ISIs. 
The multi-modal distribution consists of four unimodal dis- 
tributions characterized by two successive firing modes; for 
example, S-D denotes the distribution of the ISI between the 
initial singlet and the final doublet firings, (c) Comparisons 
between the WSTA-estimated PRCs and the theoretical ones, 
(d) A distinct difference between two WSTA-estimated PRCs 
leads to a transition in the synchronization for a two-neuron 
system. 



randomly, under the influence of the injected fluctuating 
current, as shown in FigHJa). For the bursting neuron, 
here we adopt the chattering neuron model, which ex- 
hibits an increasing number of intra-burst spikes, such as 
from singlet to doublet, when the injected current is in- 
creased |9(. In the case considered in Fig|4] the injected 
current fluctuates about the mean level corresponding to 
the transition between singlet and doublet firings. Using 
the WSTA with an additional procedure, we can simul- 
taneously reconstruct two PRCs corresponding to singlet 
and doublet states from a single trial. The key idea here 
is to separately calculate two conditional WSTAs, de- 
pending on whether a singlet or doublet firing occurs. 
In other words, the PRC for the singlet (doublet) fir- 
ing can be reconstructed by calculating the WSTA over 
only the injected currents generating the singlet (dou- 
blet) firing. In FigfJJh), closer investigation reveals that 
the multi-modal distribution of the ISI can be regarded as 
a mixture of four different unimodal distributions, which 
are specified by two successive firing modes. From each 
unimodal distribution, the average ISI is separately com- 
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FIG. 5: The optimal strengths of the injected fluctuation in 
a noisy environment. As the background noise, zero-mean 
O-U noise with a = 0.001 was added to the dynamics of 
the potassium activation variable. In all cases, the WSTAs 
were calculated over 1000 spikes, (a-b) The estimated PRCs 
for two types of Morris-Lecar models: (a) type I excitability 
near the saddle-node bifurcation point (/bias = 8.6^iA/cm 2 , 
T- 1 =28Hz); (b) type II excitability near the Hopf bifurca- 
tion point (/bias = 30/iA/cm 2 , T _1 =64Hz). For smoothing 
algorithms, we considered several statistical regression mod- 
els, in which Fourier (F) and von Mises (M) functions were 
used as basis functions [llfl . (c) The dependence of the root- 
mean-square error (RMSE) on cr/ibias for the type I model (5 
trials), (d) The same graph as in (c) for the type II model. 



puted for calculation of the A< in Eq.(4). Finally, two 
rescaled sets of data with the same final firing mode are 
used to obtain the two PRCs for the singlet and dou- 
blet modes. Figure 0|c) demonstrates that both of the 
PRCs obtained in this way are in reasonably good agree- 
ment with the exact ones. We emphasize that two differ- 
ent PRCs can be reconstructed from a single-trial data 
set, even when the naive conventional method employ- 
ing a pulse-like perturbation is not practical. Noting 
the sharp difference between these two WSTA-estimated 
PRCs, we predict a transition between anti-synchronous 
and synchronous firing for a two-neuron system with ex- 
citatory synaptic couplings, as shown in a previous study 
O(FigHd)). 

Discussion - In actual applications to real systems, we 
have to choose an appropriate magnitude of the fluctua- 
tion of the injected current. Since some unpredictable 
and uncontrollable noise is inevitable in real systems, 
we examine the discrepancy between the estimated PRC 
and the true one in the presence of unknown background 
noise, as summarized in Fig. [5J We find that for two 
types of Morris-Lecar models near the threshold for fir- 
ing , the WSTA estimates the PRC most precisely in 
the case that an intermediate magnitude of the fluctua- 
tion of the injected current is chosen. This is because the 
signal of the PRC becomes lost in the background noise 
if the injected fluctuations are too weak, while if the fluc- 
tuations are too strong, the approximation of the phase 
description becomes poor. Figures HJc) and (d) suggest 
that the same results hold for all types of bifurcations 



and smoothing algorithms considered. 

We now give some comments on three relevant studies. 
First, when the fluctuation of the ISIs vanishes, all the cl 
weights of the WSTA are the same. This is equivalent 
to the situation in which the STA yields Z'(<p) instead of 
•^(0) 0- This superficial inconsistency can be resolved 
by considering the fact that the difference between the 
WSTA and the STA converges to zero as /x 2 . Second, cl 
recent studies have pointed out that a correction term 
appears in the phase description ([3]) if the limit-cycle 
oscillation is perturbed by the noise with a correlation 
time shorter than the relaxation time of the limit-cycle 
attractor[l3j]. Under real conditions, a fluctuating noisy c2 
signal can be regarded as a smooth signal in the limit 
of a small time scale. Therefore, Eq. ([3]) and our re- 



sult are practically valid. Third, in another study [l4|. cl 
the collective PRC to a macroscopic external force for 
globally coupld oscillators is investigated. In principle, c2 
our method can be applied to measure such a collective 
PRC. Closer examinations of the above-mentioned issues c3 
are beyond the scope of this paper, but they are of great 
importance and should be studied further. 

In conclusion, we demonstrated that the PRC can be 
reconstructed from our proposed WSTA, in which the 
weight is proportional to the magnitude of the fluctua- 
tion of the period. In this study, considering limit-cycle 
oscillations observed ubiquitously in nonlinear dissipa- 
tive systems far from equilibrium, we found a theoreti- 
cal relation between the fluctuations in the system and 
its response to an external force. This might provide 
useful insight for further studies as the development of 
fluctuation dissipation theorem. In comparison with the 
standard method employing pulse- like perturbations, fur- 
thermore, the proposed WSTA method is experimentally 
reliable, fast and wide applicable, as shown by the fact 
that two different PRCs can be reproduced even from a 
single-trial data set for two-mode bursting neurons. We 
believe that this method will contribute to elucidating 
the nature of real dynamical system in a broad range of 
contexts. 
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